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ABSTRACT 

In the sequential accretion model, planets form through the sedimentation of dust, cohesive 
collisions of planetesimals, and coagulation of protoplanetary embryos prior to the onset of ef- 
ficient gas accretion. As progenitors of terrestrial planets and the cores of gas giant planets, 
embryos have comparable masses and are separated by the full width of their feeding zones after 
the oligarchic growth. In this context, we investigate the orbit-crossing time (Tc) of protoplanet 
systems both with and without a gas-disk background. The protoplanets are initially with equal 
masses and separation (EMS systems) scaled by their mutual Hill's radii. In a gas-free envi- 
ronment, log(Tc/yr) ~ A -I- i? log(fco/2.3), where fco is the initial separation of the protoplanets 
normalized by their Hill's radii, A and B arc functions of their masses and initial eccentrici- 
ties. Through a simple analytical approach, we demonstrate that the evolution of the velocity 
dispersion in an EMS system follows a random walk. The stochastic nature of random-walk 
diffusion leads to (i) an increasing average eccentricity < e >(x t^/^, where t is the time; (ii) 
Rayleigh-distributed eccentricities {P{e,t) = e/cr^ exp(— e^/(2a^)), where P is the probability 
and a{t) is the dispersion) of the protoplanets; (iii) a power-law dependence of Tc on planetary 
separation. As evidence for the chaotic diffusion, the observed eccentricities of known extra solar 
planets can be approximated by a Rayleigh distribution. In a gaseous environment, eccentricities 
of the protoplanetary embryos are damped by their interactions with the gas disk on a time scale 
Ttidai which is inversely proportional to the surface density of the gas. When they become well 
separated (with fco — 6 — 12), the orbit-crossing tendency of embryos is suppressed by the tidal 
drag and their growth is stalled along with low-eccentricity orbits. However, the efficiency of 
tidal damping declines with the gas depletion. We evaluate the isolation masses of the embryos, 
which determine the probability of gas giant formation, as a function of the dust and gas surface 
densities. Similar processes regulate the early evolution of multiple gas giant planet systems. 
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1. Introduction 

The origin and evolution of multiple planets 
around the Sun have been the primary stimu- 
lus for the studies of classical N-body systems. 
The pioneering analysis of Poincare (1892) estab- 
lished the paradigm that the motion of systems 
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with more than two bodies is not integrable. This 
fundamental result is well supported by the dy- 
namical diversity among the ~ 200 extra solar 
planetary systems discovered in the past decade. 
While considerable efforts have been made to in- 
terpret various physical processes which may con- 
tribute to their new-found properties, it may be 
fruitful to explore, in depth, the implication of 
long-term dynamical interactions among members 
of multiple-planet systems. These analysis may 
eventually provide the basis for a theory of statis- 
tical mechanics which characterizes the architec- 
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ture of planetary systems. 

Statistical mechanics has been employed to 
study other N-body systems in astrophysics. In 
the context of stellar clusters, the time scale of 
phase-space relaxation may be evaluated by a 
Fokker- Planck approximation. The magnitude 
of the diffusion coefficient is determined by an 
impulse approximation, i.e. as an ensemble of 
independent close encounters. But in planetary 
systems, the host stars dominate the gravity field. 
Although planetary perturbations are weak, they 
persist and are correlated over many orbits. This 
aspect of the dynamical evolution makes the de- 
velopment of a statistical approach particularly 
difficult. 

The investigation of the phase space diffu- 
sion is closely related to the stability of plan- 
etary systems. With the exception of periodic 
and quasi-periodic orbits, the stability of most 
orbits in a general N-body planetary system 
is not known. The Kolmogorov-Arnold-Moser 
(KAM) theory proved that a non-degenerate 
integrable Hamiltonian system may preserve 
most of its stable (quasi-periodic) motions un- 
der sufficiently sma ll and analytic al perturbation s 
(jKolmogorov 1954 iMoser IQSSt lArnold 1963^. 
For those non-stable motions, the Nekhoroshcv 
theorem showed that, the time that an orbit 
becomes unstable grows exponentially with re- 
spect to the inverse of th e non-integrable param- 
eter (jNekhoroshev 19771 ). For vanishing "per- 
turbation" amplitude, the diffusion time scale 
become infinitely long. However, most systems 
of astronomical interest, such as planetary sys- 
tems, are degenerate. Consequently, the appli- 
cations of the powerful KAM and Nekhoroshev 
theorems turned out to be indirect and difficu lt 
(|Siegel fc Moser 197ll : iMorbideUi fc Cxuzzo 19971 ). 

Nevertheless, the stability of planetary sys- 
tems remains an important problem with many 
applications. The first application of this fun- 
damental issue concerns the dynamical age of 
the Solar System. Although interactions be- 
tween the planets give rise to chaotic motions, 
the system is expected to remain essentially sta- 
ble over a time much longer than its present age of 
4.6 Gvr dLask ar 1989^: [Sussman fc Wisdom 199^ 
Murrav fc Holman 1999 1 

Another issue is the stability of a proto-planet 
system during the early stage of its formation. Ac- 



cording to the conventional sequential-accretion 
scenario, the terrestrial planets are formed by the 
coagulation o f planctesi mals in protostellar disks 
(ISafronov 19691 : IWetherill 1980l) . Through sev- 
eral stages of runaway and oligarchic growth, co- 
hesive collisions lead to the emerge nce of mas- 
sive protoplanet ary embryos (fKokubo fc Ida 20021 : 
Ida fc Lin 20041). According to the numerical sim- 



ulations (Kokub ofc Ida 1998h . protoplanets form 
with comparable masses and similar separation 
(~ 10 Hill's radii). The stability of such proto- 
planet systems could be crucial for the subsequent 
evolutions and final configurations of the system, 
like the presence of Earth-mass planets near their 
host stars (e.g., Zhou et al. 2005). 

A third issue concerns the excitation of the 
large eccentricities as well as the stability of the re- 
cently observed extra solar planet system^ The 
observed extra solar planet systems have a median 
eccentricity of 0.25 ( Marcv et al. 2005h . Despite 
its large uncertainties, the eccentricity distribu- 
tion of extra solar planets is quite different from 
our Solar System. As interactions between gaseous 
disks and protoplanets a re expected to generall y 
limit their eccentricities ( Papaloizou et al. 2006f ). 
the origin of the large eccentricities in extra solar 
systems remains poorly understood. 

Despite these important questions, an analytic 
theory for stability of planetary systems has not 
been attained. Facing this enormous complex- 
ity, recent attempts to understand some aspects 
of this process have been reduced to a subset of 
three-body problems. Based on the results from 
qualitative studies of the general three-body prob- 
lem (e.g., Marchal 1990), Gladman (1993) inves- 
tigated the stability of the two planet systems 
both analytically and numerically. He found that 
a system of two planets with mass ratios to the 
star /ii , /i2 could be Hill stable if their separation 
> 2A/3(^i4^)^^^' ^'^'^^^ HiU stable is defined as 
orbits that will never cross. In systems with more 
than two planets, the most practical approach is 
to resort to numerical simulations. Due to the 
large degrees of freedom of these systems, restric- 
tions are needed to reduce the range of configu- 
rations for parameter studies. Motivated by the 
characteristics of embryo systems after runaway 
and oligarchic growth, a series of investigations 
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have been carried out to study idealized but well- 
defined planetary systems with equal masses and 
scaled separation. Hereafter we refer these ideal- 
ized planet systems as EMS systems. 

Chambers et al. (1996) determined numerically 
the orbital crossing time Tc of EMS systems with 
n planets (n > 3) initially on circular orbits. Tlicy 
found an exponential relation logTc ~ ko, which 
seems to be independent of n. The dimension- 
less parameter kg is the scaled initial separation. 
They did not provide any explanation of the un- 
derlining cause of this relation. Later, Yoshinaga, 
Kokubo and Makino (1999) generalized this study 
to the cases that the planets are initially on non- 
circular and non-coplanar orbits. In the limit of 
small initial eccentricity eo and inclination, they 
obtained similar results as previous investigators. 
Later, the instability of EMS systems under so- 
lar nebular gas drag was studied by Iwasaki et al. 
(2001, 2002) and Iwasaki & Ohtsuki (2006). 

However, the EMS systems studied in these 
works are with separation ko < 10. For realistic 
planetary systems, the initial separation between 
planets may be larger, with a gas disk during the 
stage of planet formation. In the Solar System, 
the present-day values of fco ^ 8 — 64. According 
to the numerical simulations of planet formation 
(Kokubo & Ida 2002, Ida & Lin 2004), after the 
planetary embryos have depleted nearby planetes- 
imals and reached isolation masses, the embryos 
were separated with fco ^ 10 — 12. 

The initial motivation of the present work is to 
extend the previous studies to the cases fco > 10 
both with and without a gas disk, and to de- 
rive a functional dependence of Tc on fco,/i, eg. 
We show in §2 that, the orbit crossing time Tc 
is better approximated by a power-law relation 
log Tc ^ log fco . A simple analytical interpreta- 
tion of this relation is suggested in §3. We also 
show that the average eccentricity of an EMS sys- 
tem in a gas- free environment increases as ^ t^^"^ . 
We identify this evolution as a result of the ran- 
dom walk diffusion in phase space which accounts 
for the power-law dependence of the orbital cross- 
ing time on the initial separation. In §4, we ex- 
tend the study to the cases when the protoplan- 
ets (or embryos) arc embedded in a gas environ- 
ment. This investigation determines the range of 
feeding zones and isolation masses of embryos in 
gas-rich protostellar disks. The embryos' masses 



and separations during the post-oligarchic evolu- 
tion in a depleting gas environment are derived. 
These quantities determine the probability of gas 
giant formation. We show that the observed ec- 
centricity distribution of known extra solar plan- 
ets has the form of a Rayleigh distribution. We 
cite this property as evidence for chaotic diffusion 
being the dominant excitation mechanism. Sum- 
mary and the implications of our results on the 
formation of planet systems are presented in the 
final section. 

2. Empirical formula for Tc without gas 
disk 

The model of an EMS system is given as follows. 
Suppose n protoplanets (or planets for simplicity) 
with equal masses move aroimd a star with one 
solar mass, and the separation between them are 
equal when scaled by their mutual Hill's radii. In 
this paper all the orbits of the planets are copla- 
nar, especially the EMS systems are in a gas-free 
environment in this and the coming sections. 

We denote the mass ratios of the planets to the 
star, the semi-major axes and eccentricities of the 
planets' orbits as /x, aj and (i=l,...,n), respec- 
tively. The scaled separation and eccentricities of 
the planet orbits are 

^ ~ ' iijr ' ' (» = 1> •••>"■- 1), Q\ 

el = f, ii = 1, ...,n), 

respectively, where Rh is the mutual Hill's radius 
and h is the relative separation of two neighboring 
planets, defined as 

= (2M^i/3«i + «i+i^ ^ ^ "m-"-'. (2) 
3 2 tti+i + Qi 

Thus the orbits of two neighboring planets with 
e = 1 will cross if the difference between their 
perihelion angles is tt. For simplicity, we adopt 
the same initial eccentricities cq, while the initial 
mean anomaly Mi, (i = 1, ...,n), and longitude of 
perihelion vui of each planet are chosen randomly. 
We take n = 9, and arbitrarily specify the initial 
semi-major axis of the fourth planet a4 = lAU 
for normalization purposes. So when the initial 
separation ko = k{t = 0) varies, the planet system 
is enlarged both inward and outward. 

The orbital crossing time of the EMS system 
(denoted as Tc) is defined as the minimum dura- 
tion when either of the following two situations 
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occurs between any two planets during the evolu- 
tion: (1) close encounter, defined as the distance 
between them is less than their mutual Hill's ra- 
dius, (2) orbit crossing, defined as Ui > a^+i, (i = 
1, ...,n— 1). We use the symplectic code of Wis- 
dom and Holman (1991) from the SWIFT package 
(Levison & Duncan 1994). Whenever orbit cross- 
ing or a close encounter occurs, we halt the inte- 
gration. The time step is chosen to accommodate 
~ 20 steps per inner planet orbit, and the accu- 
mulated error of the relative energy during the in- 
tegration is constrained to be ^ 10~^° — 10~^ until 
the system becomes unstable. 

We investigate mainly 7 typical values oi fj, = 
10*, {i = — 10, —4). For each value of /i, we do 
10 sets of simulations with initial eccentricities of 
the planets in the range e = 0, 0.1, 0.2, 0.9. For 
each set of parameters, many orbits with various 
initial value fco are integrated to determine the re- 
lationship between Tc and fco- 

Fig.l shows the dependence of Tc on fco for a 
range of fi. We find there exists roughly a critical 
kc such that, Tc is independent of fco for fco < 
fcc and increases with fco for fco > fcc(Fig.la,lb). 
These two branches of solutions join continuously 
at fco = fcc with the approximation Tc{ko ~ fcc) = 

A. We are primarily interested in the range of 
fco > fcc for which the numerical results can be 
fitted with log(Tc/yr) = A-l-B log(ko/kc). In order 
to obtain the value of the numerical coefficients. A, 

B, and fcc, we proceed as follows: 

(i) We first determine fcc by scaling Tc with 

fco in the range [1.5,3.5]. We found the 
eccentricity-dependence of fcc to be negligi- 
ble over e G [0,0.5]. For the entire range of 
/i, we obtain fcc ~ 2.3, again insensitive to 
the magnitude of /i (Fig. 2a). 

(ii) We evaluate the average values oi A = 

Tc(fco = fcc), and find A = (-0.91 ± 
0.08) - (0.27 ± 0.01)logAf (Fig.2b). A 
more general expression, which also incor- 
porates the eccentricity dependence of T^, is 
A = -2 + go - 0.27 log /i. 

(ill) Finally, we determine the magnitude of B. 
From the slopes of the log(rc) — log(fco) 
curves of Fig.l, we obtain the eccentricity 
and /i dependence of B (Fig.2c-d). A rea- 
sonable approximation for the B{fi,eo) is 



B = bi + 62 log M + {h + h^ogiJ,)eo, with 
61 = 18. 7 ±0.6, 62 = 1.11 ± 0.08, 63 = 
-16.8 ±0.6, 64 = -1.24 ±0.08. 

After some exhaustive simulations, we obtain the 
following empirical fitting formula: 

log(§) = A + i?log(i%). 

(fco > 2.3,10-4 ^Q-io) y-^) 

where 

A= (-2 + go -0.27 log ^) 

B = (18.7+ 1.1 log /i) - (16.8 + 1.2 log Ai)eo. 

(4) 

The predictions given by the formula ([3]) are 
plotted also in Fig. 1. Wc find the formula 
agrees well with the numerical results for plan- 
etary masses lO""' < jJ. < 10"^*^. In this mass 
range, slope B is positive. The above formula 
(|3l) generalizes a similar approach introduced by 
Chambers et al. (1996jl. The distribution of Tc 
in the separation-mass (fco — /i) space is shown in 
Fig. 3a for eo = 0. 

However, we find formula ([3]) is not satisfied 
when applied to ^ ~ lO-''. Since in these sit- 
uations, resonances between planets are strong 
and dominate the dynamics at the place fco = 
2(f^T)/(|M)^/^ where q = {n./n.+ifl^ is the ra- 
tio of the mean motions of planets i and z + 1 . As 

~ 10~^ is the ideal case for giant planet sys- 
tems, we investigate this case for planets on initial 
circular orbits, and find the orbital crossing time 
can be approximated by a simple formula in the 
case fco < 10: 

log(^) w -5.0 + 2.2fco. (m 10"^ g = 0) (5) 

Fig. 3b shows the numerically determined orbital 
crossing time with the best fit formula The 
drop of Tc near fco ~ 5 is due to the presence of 
the 2 : 1 resonance (fco — 5.2) between the planets. 



^For eo = and fi = 10~^, Chambers et al. (1996) found 
logTc = fefco + c in the range fco < 10, with b = 0.76 ± 
0.03 and c = —0.36 ± 0.18. They also obtained similar 
expressions for other values of /i. This expression can be 
obtained from equation (|3} in the limit of small fco. For 
example, in the range of fc < 10, x = (fco — 6)/6 < 1 and 
equation l(3)l reduces to logTc = ll[log(l + x) + log(2^)] — 
0.11 « -r^x + 4.47 = O.SOfco - 0.31. 
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Fig. 1. — Variations of the orbit- crossing time Tc with 
initial orbital separation fco in the 9-planet EMS sys- 
tems of different ^ and eo. The triangles, squares and 
crosses denote systems with eo = 0, 0.5, 0.9, respec- 
tively. The solid lines are calculated from the empiri- 
cal formula ([Sjl. In the ^ = 10~® case (d), a correction 
of +0.5 is added to the values of log Tc given by equa- 
tion Q. 




Fig. 2. — The procedure to determine the coefficients 
kc,A,B in formula ([3]). (a) Variations of the average 
Tc with small fco. The average is taken over e G [0, 0.5]. 
From bottom to up, the curves correspond to EMS 
systems with fi = 10^*, 10~', respectively, kc is de- 
fined so that < Tc > begins to increase with fco at 
fco > fcc (b) Determine A =< Tc > {k = kc) for 
different fi. The squares with error bars are numer- 
ical results, while the solid line {A — Ai + A2 log fi) 
is the best-fit line. The best-fit coefficients are also 
shown, (c) The triangles, squares and circles with 
error bars denote the best-fit slopes B of the curves 
(log(Tc) — log(fco)) in Fig.l. As a function of eo, it 
can be expressed as _B = Bi + -6260 for various fi. 
The best-fit coefficients for Bi — bi + 62 log(/i) and 
B2 = bs + &4 log(At) are shown in (d). 



From equation ([3]), we can highlight the differ- 
ence in the crossing time of two EMS systems (de- 
noted as SI and S2, respectively) on initial circular 
orbits: 

• Suppose SI and S2 have the same planetary 
masses: ^1 = 112 = 



^cl _ ^ ^01 -|18. 7+1.1 loRfj. 
Tc2 ko2 



(6) 



Thus for example, if = — 7 and fco 1/^02 = 
2, the above formula yields Tci/Tc2 « 2000. 
The crossing time of the widely separated 
system (SI) is three orders of magnitude 
larger than that of the compact system (S2), 
even though the initial separation among 
planets differs only by a factor of 2. 

• In contrast, let SI and S2 have the same 
planet separation fcoi = fco2 = ^o, 



£cl _ j-M]_x-0.27+l.llog(fco/2.3) 
Tc2 m 



(7) 



Thus for example, if ko = 10 and /ii//i2 = 
10, it gives Tci/Tc2 « 2.7. The crossing time 
for the massive system (SI) is around three 
times longer than the less massive system 
(S2), provided their normalized (by the Hill's 
radius) separations are the same. 

3. A simple analytical approximation 

The numerical simulations, though informative, 
do not provide any underlying theory for the ori- 
gin of the dependence of Tc on ko, /i and cq. In 
this section, we present a simple analytical ap- 
proach in an attempt to describe the evolution of 
the EMS systems without gas disk. We identify 
the planets of an EMS system with subscript I 
(1, 2, Z — 1,1,1 + 1, n with n > 3), in the in- 
creasing order of their initial semi-major axes. We 
consider the evolution of a representative planet 
1 < I < n. Assume all the planets are initially 
on circular orbits, and in the limit of close separa- 
tion, i.e. ai+i — ai << ai. According to equations 
(HI) and this approximation is equivalent to 
fco(2/i/3)^/'^ ^ 1. We call it the close separation 
assumption. The largest perturbations on planet 
I come from close encounters with nearby planets 
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Fig. 3. — The orbital crossing time on parameter 
space, (a): Contour lines of log(rc) of EMS systems in 
circular orbits in the space of initial orbital separation 
ko and planet masses /i. The numbers in the curves 
are log(rc). They are obtained from formula (b): 
Variations of Tc on fco for /i = 10""^. Squares are from 
numerical simulations, and the solid line is from for- 
mula (O. The big drop at fco ~ 5 corresponds to 2 : 1 
resonance between planets. 
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Fig. 4. — Evolution of g = a(l — e),a,Q — a(l + e) 
for the 9-planet EMS system in a (a) gas-free, (b) gas- 
rich environment. Parameters in (a) are fi = 10~^, 
eo = 0, fco = 8. The orbital crossing time is 7 x 10^ 
yr, according to equation (|3]). Parameters in (b) are 
fj, = lO"'^, eo = 0.5/i, ko — 8. The orbital crossing 
time is 1.5 x lO^yr. From Fig. 3 and formula Q, the 
orbital crossing time for the same parameters but in a 
gas-free environment is ~ 10* yr. 



(planet l±l). Under the close separation assump- 
tion, the interactions between each pair of neigh- 
bors can be well approximated by an independent 
set of Hill's problems. 

We define e = (a; - ai-i)/ai ~ fco(2/i/3)i/3 
as the relative semi-major axis, zi = e; exp(in7/) 
as the Runge-Lenz vector, and zui as the longi- 
tude of periapse of planet I. We consider the limit 
e; <C e <C 1. To first order in /i, a;,ai_i do not 
change during close encounters (Henon & Petit 
1986). We assume that during all close encoun- 
ters prior to orbit crossing the semi-major axes of 
the planets do not have significant secular changes. 
This assumption is supported by the numerical re- 
sults (See Fig. 4a). However, z evolves and after 
the j-th close counter between the planets I — 1 
and /, the change in z is given as 

= ^i-i - exp(iAj_i), (j > 1), (8) 

where Xj is the mean longitude of planet I , g = 
|[2Jro(|) + Kiil)] ~ 2.24, where Kq and Ki are 
modified Bessel functions (Henon & Petit 1986, 
Duncan, Quinn & Tremaine 1989). The time be- 
tween two consecutive close encounters is given as 
T, = Tj[(aj/a,_i)3/2 _ l]-i « |r,e-i, where is 
the orbital period of the planet I. 

For illustrative purposes, we adopt a; = 1 AU, 
so T; = 1 yr, and the change of A during one 
encounter is given as Xj « Aj_i -I- ^ . Since 
e <C 1 and the change of e is second order in /i, 
Xj (j — l,...,rt) at successive encounters behave 
like a series of random numbers in [0, 27r]. Ac- 
cording to ([8]) we have, 

2 2 r,9M • /A ^ 5^^^ 

ej - Sj-i = -^'—e^-i sm(Aj_i - zuj^i) + 

'(9) 

Due to the near-random phase of Xj , the first term 
in equation ^ averages to zero over a long time. 
Changes of induced by the perturbations from 
planets /±2,Z±3, ... are ^ 1/2^,1/3'*,... times 
those from Z ± 1. However, the periods of close 
encounters between planet I and these planets are 
^ 1/2,1/3,..., times Tg, respectively. Therefore, 
when we take account of perturbations from more 
distant planets on both sides, we introduce a factor 
2(1 + 1/23 _^ 1/33 _|_ ) ^ 2.40, so that < Ae^ 
^Ag^ li^t'^ . The average eccentricity of the Z-th 
planet after j close encounters with nearby planets 
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is estimated to be 



yr 
(10) 

where we have substituted j = t/Ts = |ei/yr. 
This formula will be confirmed by numerical sim- 
ulations in this section. 

According to the criteria specified in §2, orbit 
crossing occurs when < >i/2^ h = ifco(|/i)^/'^. 
From equation (fTU]) . we derive, 



log(^) « -1.1 

yr 



Slogfco - ilog^. (11) 



This expression describes the power law depen- 
dence of Tc on fco as in equation However, the 
discrepancy between the coefficients B and 5 in 
equations ^ and ifTTj) is considerable, especially 
when jj, is large. This may be due to the close 
separation assumption, e ~ ko^i^^^ ^ 1 no longer 
being valid for moderate fcg and /i > 10~^. More- 
over, the sign of the coefficient of log/i is negative 
which disagrees with equation ([3]). This may be 
caused by the oversimplified assumptions in the 
analytical model. 

Next, we show that the evolution of the av- 
erage eccentricity (< >^/^(x t^^^) is mainly 
driven by a random walk process. The stochas- 
tic nature of the perturbations also leads to the 
power law dependence of Tc on fco. We define the 
velocity dispersion as w = |vkcp| — |vcir|, where 
Vkop,Vcir are the velocities of Keplerian and cir- 
cular motion respectively. It is easy to show that 
V = nae cos f + o{e^) , where / is the true anomaly. 
We consider a group of orbits in phase space, and 
the probability of planet / having velocity disper- 
sion V is denoted by P{v). Thus P{v) describes 
the distribution of a group of orbits in velocity 
dispersion space. Since every close encounter be- 
tween planets will modify the distribution, P{v) is 
a function of time t (or j encounters). We assume 
that the planetary motions are chaotic and occupy 
a stochastic region in the phase space. This as- 
sumption is justified by the random phase of A 
and the non-zero Lyapunov exponents shown at 
the end of this section. 

Under the chaotic assumption, the evolution of 
P{v,j) obeys the Fokker-Planck equation (Licht- 
enberg & Lieberman 1990): 



dP 



d_ 

dv 



{BP) 



2dv^ 



(DP), (12) 



where B, D are the frictional and diffusion coeffi- 
cients, respectively, with 



r27T 



^V^/J^^[Ae(V')cos/]2dV', 



(13) 



where ip = X — w. Following the standard pro- 
cedure in celestial mechanics, we carry out orbit 
averaging around the Keplerian motion so that 
cos^ / = 1/2-1- o(e^). We adopt the approxima- 
tion (Ae)^ « Ae^. According to equation ([9]), we 
find D{v) sa n? a? ^"^ . Since D is indepen- 
dent of w, i? = = 0. After replacing j by i, 
the Fokkcr-Planck equation is converted into the 
standard diffusion equation: 



dP _ j^d^P 
dt dv^ ' 



(14) 



where D = jeDyr ^ 



The time dependent solution of the above equa- 
tion with the initial value P(w,0) = (5(0) (where 
5{x) is the Dirac delta function) is a Gaussian (i.e., 
normal) distribution: 



P{v,t) 



Substituting _D, we find 



^cxp(-;^), a^{2btf'\ 
o\/2tt 2a^ 



na 



(15) 



(16) 



We convert equation to a distribution of ec- 
centricity by substituting v = nae cos /, where 
functions of cos / are replaced by the average val- 
ues over a Keplerian period, < cos/ >— ~e and 
< cos^ / >= 1/2. Thus we get, 



6 ^ 

P(e,t) = — exp(-— 2), a 



\/2cr 



(17) 



which has the form of a Raylcigh distribution. 

In order to verify the above analytical results, 
we carry out some numerical simulations with 
EMS systems of n = 50 protoplanets. These re- 
sults also provide a self-consistent verification on 
the assumed chaotic nature of planetary motion. 
In these simulations, we specify the following ini- 
tial conditions. The planets are initially placed 
on circular orbits, with ai = lAU. We utilize the 
Hermit scheme P(EC)^ in order to follow the plan- 
ets' evolution after their orbital crossing (Makino 
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& Aarseth 1992, Yoshinaga, Kokubo & Makino 
1999). 

Figs. 5 and 6 show some typical numerical re- 
sults. At each given epoch, the normalized veloc- 
ity dispersions relative to the circular orbits fol- 
low a Gaussian distribution (jlSp . The correspond- 
ing eccentricities obey a Rayleigh distribution (fT7|) 
(see Fig. 5). Fig. 6 shows the evolution of the nor- 
malized velocity dispersion and that of the aver- 
age eccentricity. Both quantities grow with t^/^ 
as predicted by the analytical approach in equa- 
tions (fT6|) and (fTO|) . The agreements are excellent 
for jj, = 10^^ and 10^^. Similar to the Brownian 
motion, the evolution of the velocity dispersion in 
an EMS system is a random walk process. How- 
ever, the coefficients are not well predicted by the 
analytic expression for fi = 10~^. The less satis- 
factory predictions of equations (fT6|) and pO|) for 
large masses may be due to the close separation 
assumption e ~ k^y}/^ <^ 1 being poorly satisfied 
in the limit /i > 10^'"'. We note that in Fig. 6 there 
are no very significant transitions in the evolution 
of < e > when orbit crossing occurs (~ 10"^ — 10** 
yr according to Fig. 3a). This behavior indicates 
that the growth of < e > is a result of a slow 
diffusion process. 

We now justify the assumption of stochastic 
phase space. For this task, we calculate the Lya- 
punov exponents (LE) at a finite time x(i) for 
the EMS systems. As is well established for two- 
planet systems, there is a well-defined boundary 
between the regular and chaotic motions which is 
demarcated by fco ~ 2/x2/7(Wisdom 1980, Glad- 
man 1993). However, in EMS systems with n > 3, 
xit) may undergo transitions to a finite value af- 
ter a long period of time. The reason for this be- 
havior is due to the increase of velocity dispersion 
(~ t^/^) through orbital diffusion. Orbits initially 
in a regular region will finally, though after a very 
long time, become chaotic due to the increase of 
velocity dispersion. Thus we believe the changing 
from chaotic motion to regular motion along fco 
space is gradual, and there is no clear boundary 
between the domains of regular and chaotic mo- 
tions (Fig. 7). We will discuss this problem else- 
where (Zhou & Sun 2007). In Fig. 8, we map out 
the Lyapunov time {T^, inverse of LE) as a func- 
tion of (fcoi^)- For computational simplicity, we 
consider here only those systems on circular orbits 
initially. The chaotic nature of the entire param- 




-0.5 0.0 0.5 

velocity dispersion 



eccentricity 



Fig. 5. — Distributions of (a) the velocity dispersions 

V and (b) eccentricities in four runs of 50-planet EMS 
systems with fj, = 10*, fco = 5 at time t = 0.4 Myr. 
The fit Gaussian distribution in (a) is according to 
equation (|15p with a = 0.336, an adjustment of < 

V >= -0.0342, and a scale factor of 37.4. The fit of the 
Rayleigh distribution in (b) is according to equation 
(fT7)l with cr = 0.194 and a scale factor of 10. 



(a) 





















(b) 


A 








/ ^^"^^ c 



log t (yr) 



log t (yr) 



Fig. 6. — Evolution of (a) the variances of velocity 
dispersions a normalized by na and (b) the average 
eccentricities in a 50-planet EMS system with fco = 
5 and different ^i: A. = 10"^ B. ^ = 10"^ C. 
/i = 10~^. n, a are the mean motion and semi-major 
axis of each planet. The solid lines in (a) and (b) are 
obtained from the analytical formulas (|16|l and (I10|l . 
respectively. 
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eter domain calculated justifies our random-phase 
assumption. 

We also plot in Fig. 8 three lines of constant Tc 
derived from equation ([3]). The line corresponds 
to Tc = 10^'^ yr lies on the boundary between 
the strongly (with Tl < 10'^ yr) and weakly (with 
Tl > 10^ yr) chaotic regions. In comparison with 
Fig. 4, we find, that the Luapunov time of an EMS 
system in the strongly chaotic region is essentially 
independent of fco, while in the weakly chaotic re- 
gions, Tl is correlated with Tc, large Tc implies 
large T^. This indicates that the Lyapunov time 
can be either correlated with or independent of the 
orbital crossing time, which is a counter example 
to the conjecture proposed by Lecar et al. (1992). 

4. Presence of gas disk 

As indicated in the abstract and introduction, 
one motivation for our present study is to consider 
the growth of protoplanetary embryos as they un- 
dergo a transition from dynamical isolation to 
post-oligarchic evolution. The above analysis on 
the evolution of EMS systems in a gas-free envi- 
ronment is appropriate for late stages after the gas 
depletion. In this section, we consider the stability 
of EMS systems in a gas environment. Intuitively, 
gas provides an eccentricity damping mechanism 
which may suppress the growth of velocity disper- 
sion and thus prolong the orbit crossing time. 

For illustration, we adopt a fiducial model for 
the gas surface density based on the minimum 
mass nebula model such that 



- So/g/dcp(Y^) 



-3/2 



(18) 



where S 



2400gcm ^ and is a sc aling factor 



(|Havashi et al. 19851 : llda fc Lin 200l . We also 
use an idealized prescription to approximate the 
decline of the gas surface density with a uniform 
depletion faction /dcp = exp(— i/Tdcp). We adopt 
a magnitude for the gas depletion time scale to be 
Tdep = 3 Myr based on observations (Haisch et al. 
2001). 

In a gaseous disk background, a protoplanet 
with mass ratio fi suffers a gravitational tidal drag, 
which for simplicity, can be expressed as 



Ftidal - -T^^idA^ - ^^c), 



(19) 



where V and Vj. are the Keplerian and circular ve- 
locity of the protoplanet, respectively (Kominami 




2.5 3.0 3.5 4.D 4.5 5.0 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

logt (yr) log t (yr) 



Fig. 7. — Lyapunov exponents for orbits with fco = 
2.0 -I- i * 0.3, i = 0, 19 and fi = 10~^, eo = in an 
EMS system with (a) 2 planets, (b) 9 planets. The 
Lyapunov exponents are calculated from the varia- 
tional equations along the solutions. There are 20 lines 
in each plot which correspond to i=0,...,19. The accu- 
mulated value of relative energy error is ~ 10" ^° for 
the simulations. 




Log (TL/yr) 

u 

Fig. 8. — Lyapunov time, log(rz,), in the parame- 
ter space (fco,log(/^)) of 9-planet EMS systems with 
eo = 0. The three dashed lines A,B,C correspond to 
the crossing time of 10", 10''■^ 10^ yr, obtained from 
equation Q, respectively. 
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& Ida 2002, Nagasawa et al. 2005). The time scale 
Ttidai is defined as (Ward 1989, Artymowicz 1993) 

Ttidai « 0.75 X lO-3f-^f-^l^-\-^f yr. (20) 

For example, the magnitude of Ttidai is ^ 10** yr 
for a protoplanet with mass ratio ^ = 10^^. 

In principle, an imbalance between the tidal 
force on either side of the protoplanet 's orbit can 
lead to "type I" migration (Goldreich & Tremaine 
1980, Ward 1997). But the efficiency of this pro- 
cess may be suppressed by turbulence and non- 
linear response in the disks (KoUer et al. 2003; 
Laughlin et al. 2004; Nelson & Papaloizou 2004). 
We neglect the effect of type I migration. However, 
under the tidal force, eccentricity and inclination 
damping can also lead to semi-major axes evolu- 
tion. To the leading orders of e and i we have, 

h<w> = -«7^(5e2 + 2i2), 

=l<ll>^ L_ (21) 

e dt i dt Ttidai ' 

The relative importance of eccentricity excita- 
tion by planetary perturbations versus tidal damp- 
ing can be estimated by comparing Tc with Ttidai- 
As the damping process proceeds in an exponen- 
tial fashion, the growth of eccentricity is through 
diffusion, which does not have a distinct charac- 
teristic time scale itself. However, it has a relevant 
time scale of Tc when orbital crossing is reached. 
In addition, Ttidai oc S^^- During gas depletion, 
Ttidai increases as /dep vanishes and the efficiency 
of tidal damping weakens. On general grounds, we 
anticipate several possible limiting outcomes: 

(i) For closely-separated protoplanets, planetary 

perturbations are more effective than tidal 
damping, so we expect Tc Ttidai, and or- 
bital crossing occurring before the disk is de- 
pleted. 

(ii) In the range of modest separation, the pro- 

toplanets' eccentricities excited by their mu- 
tual interactions are effectively damped by 
the disk gas. Orbital crossing occurs only af- 
ter severe gas depletion such that Tc > Tdcp. 

(iii) Due to its very long excitation time scale 
even without a gas background, the eccen- 
tricities of widely separated protoplanets 
cannot be excited before the gas is severely 
depleted. Thus Tc is unaffected by the tidal 
damping. 



In order to verify these conjectures, we carry 
out a new set of numerical calculations, taking into 
account the tidal dissipation effect. We adopt a 
representative value /x = 10~^. In Fig. 9, we com- 
pare the results of these calculations with those 
obtained for EMS systems without any gas. 

In systems with eg = and fco < 5, Tc is not 
affected by the presence of the disk gas. Accord- 
ing to the above classification, we consider these 
systems as closely separated. However, the pres- 
ence of gas disk delays the crossing time of planets 
with modest separation (e.g., b < < S m the 
case of eo = 0) until gas depletion. Widely sep- 
arated systems (with fco > 8) are not affected by 
the presence of the gas. 

To illustrate the dominant effect of tidal drag, 
we study the evolution of an EMS system during 
the depletion of the gas disk. In Fig. 4b, we plot 
the evolutions of periapse distance q — a(l — e), 
semi-major axis a, apoapse distance Q = a(l -I- e) 
of an EMS system with modest separation (fco = 8 
and eo — 0.5). Evidently, the eccentricity growth 
occurs only after gas depletion for this system. Al- 
though the magnitude of Tc ^ 10** yr in a gas-free 
environment (Fig. 9 and eq. [5]), the tidal damp- 
ing effect prolongs it to ~ 10^ yr. 

During the epoch of oligarchic growth, embryos 
have similar masses 

~ 27rl]d(ai+i - ai)ai/M,, (22) 

where Sd is the surface density of the planetesi- 
mals and is the stellar mass. From equations 
([T]) and ([2]), we obtain 

For illustration, we adopt the surface density of a 
planetesimal disk as 

Sd - IW.U^r^'h cm"2, (24) 

where /d is a scaling constant relative to that of 
the minimum mass nebula, /ice is the volatile ice 
enhancement factor (/ice = 1 for a < 2.7 AU and 
/ice = 4.2 for a > 2.7 AU). Substituting it into 
equation (p3|) . we obtain the isolation mass, which 
depends on fco- 

Afi,o = 0.51 X 10-2Me?7fco/', (25) 
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where 



V = (/d/icc)'/'( 



'3/2 



(26) 



During the formation of protoplanets, orbital 
crossing induces protoplanets to undergo cohesive 
collisions, mass growth, and increasing separation. 
This stage corresponds to case (i). Prior to the gas 
depletion, the value of fco for an EMS system in- 
creases until the perturbation between protoplan- 
ets can no longer dominate their tidal interaction 
with the disk. During this end stage, which cor- 
responds to case (ii), the evolution of e, ^, and fco 
becomes stalled in a gas-rich environment. Until 
the gas is severely depleted, the embryos attain 
an isolation mass, which can be derived from the 
condition that ~ T^cp. Substituting this condi- 
tion with Tc from equation ([3]) for circular orbits 
(e = 0), and using the isolation mass determined 
from equation psp . we get the critical separation 
of an isolation mass: 



log(fciso) - yJh^ + OMc - b, 



where 



2.8 
3.6 



0.33 log 77, 
0.67 log 77 - 



logTdop, 



(27) 



(28) 



and T] is defined in equation (^5]) . In Fig. 10, 
we plot fciso as a function of /d and Tdop at lAU 
around a solar-type star. These results indicate 
that fciso decreases slightly with the increase of disk 
mass, which is consistent qualitatively with the 
numerical results of Kokubo and Ida (2002) . The 
isolation separation fciso and isolation mass Miso 
of the planets are plotted in the whole disk region 
for different T^ep (Fig. 11) and U (Fig. 12). For 
Tdop — 3 X 10^ yr and /d = 1, the isolation mass of 
embryos is ~ O.13M0 and their critical separation 
fciso — 8.7. These results support the assumption 
that isolated embryos are separated by a distance 
that is approxim ately ten times their Hill's radii 
(|lda fc Lin 2004 . 



5. Conclusions and applications 

In this paper, we extend the study on the or- 
bital crossing time (Tc) of n-planet systems with 
equal planetary masses and separation (EMS sys- 
tems), which was investigated by Chambers et al. 
(1996) and Yoshinaga et al. (1999). We find T^ 




6^=0, con 
e^=0.5h, con 
e^=0.9h, con 
- e^=0, dis 
-e^=0.5h, dis 
-e =0.9h, dis 



0.2 0.4 0.6 O.i 



1.0 1.2 
log k 



1.4 1.6 1.8 2.0 2.2 



Fig. 9. — Variations of the orbit-crossing time Tc with 
initial orbital separation ko in the 9-planet EMS sys- 
tems with a gas-free environment (dots, denoted by 
'con') or a gas-rich environment (curves, denoted by 
'dis'). Three sets of initial eccentricities are plotted 
for both cases, h is the relative separation defined in 
equation ([2]). 




■ (b) 


a=1AU 






f^=0.1 ... 






fj=10 


6 


7 



logTdep (yr) 

Fig. 10. — Variations of isolation separation fciso (in 
the unit of Hill's radius, defined in eq. [2]) with (a) 
disk enhancement factor and (b) gas depletion time 
scale Tdc-p at \AU . fciso is calculated from equation 
(|27[) . At Tdep = 3Myr and /d = 1, which corresponds 
to a surface density lOg cm~^ of dust at lAU, the iso- 
lation separation « 8.7 Hill's Radius and the isolation 
mass « O.lSMm. 



(a) 




=1Myr 




T, 


=3Myr 




T, 


=10Myr 


fd=1 







(AU) 




(AU) 



Fig. 11. — Variations of (a) isolation separation fciso 
and (b) isolation masses Miso with radial distance to 
the star for disk enhancement factor /d = 1 and dif- 
ferent gas depletion time scale Tjep. fciso and Miso are 
calculated from equations ([27t and (f25)) . respectively. 
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of EMS systems can be formulated as a power law 
in equation The results have the following 

implications: 

(i) The onset of instability in an EMS sys- 
tem mainly depends on the initial separation (ko). 
A qualitative inspection of equation ([3]) indicates 
that doubling /cq can enlarge T^. by several or- 
ders of magnitude. In two systems with identical 
ko, Tc increases with the planetary masses. This 
counter-intuitive result is due to the mass depen- 
dence of the planetary Hill's radii. For constant ko 
values, the un- normalized physical separation be- 
tween planets, i.e. at^i — ai, increases with their 
masses. 

ii) In a protostellar disk, a large population 
of low mass planetesimals emerge quickly. Dur- 
ing the early stage of disk evolution, the crossing 
time of planetesimals is relatively short. So the 
planetesimals will collide, merge and grow, lead- 
ing to the decline of their number density. Equa- 
tion p3)) suggests that ko of embryos increases 
with n. Since Tc increases rapidly with ko, the 
eccentricity growth due to dynamical diffusion is 
slowed down. In a gas-rich environment, the ec- 
centricities of embryos are also damped by their 
interaction with the disk gas. With mass distri- 
bution comparable to that of the minimum mass 
nebula, tidal damping becomes effective when em- 
bryos merge into bodies separated by /cq > 5. As 
the orbits of embryos are circularized, their growth 
is stalled. This result is supported by the simula- 
tions of planetesimal growth in a minimum mass 
environment, which leads to embryos with asymp- 
totic masses of ^ 10^^ g on nearly circular orbits 
with separation ~ 1 times of their Hill's radii 
(iKokubo fc Ida 19981) . 

iii) The gas accretion rate from protostellar 
disks onto their central stars decreases exponen- 
tially on a cha racteristic time scale of ~ 3 x 10^ yr 
(jHartmann 19 98). Presumably the magnitude of 



Ttg also decreases on a similar time scale, hence the 
tidal damping would become less effective. Sub- 
sequently, dynamical equilibria (in which Tc ~ 
Ttidai) are maintained with increasing separation, 
ko, while embryos merge, grow, and space out, al- 
beit at a much slower pace. When the disk gas 
is severely depleted within a few depletion time 
scales, Ttidai becomes large compared with Tdop 
and the embryo-disk interaction is no longer effec- 
tive. In a disk with minimum mass nebula (/d = 




(b) 








f^=10 








-'f^=0.1 






T*p-3 Myr 



(AU) 



a (AU) 



Fig. 12. — Variations of (a) isolation separation 
fciso and (b) isolation masses Miso with radial dis- 
tance to the star for different disk enhancement 
factor /d. fciso and Afiso are calculated from equa- 
tions ([27)1 and ((25)) . respectively, fd is disk en- 
hancement factor and Tdep = 3Myr is the time 
scale of gas depletion . 
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N.139,e>0.05 



eccentricity 



eccentricity 



Fig. 13. — Eccentricity distribution of the 139 ob- 
served extra solar planets with eccentricities > 0.05 
(from the data of Butler et al. 2006). The average 
eccentricity of these 139 planets is < e >= 0.31. (a) 
The histogram of the distribution in eccentricity. The 
solid line is the fit of a Rayleigh distribution by equa- 
tion (|17p with a — 0.25 and a scaling factor of 12.6. 
(b) The corresponding accumulative distributions for 
the observed 139 planets with e > 0.5 (dotted line) 
and for the best-fit Rayleigh distribution (solid line). 
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1), the isolation separation (fciso) and isolation 
mass (Miso) of embryo determined by Tc ~ Tdop 
are 8.7 Rh and 0.13 at 1 AU, respectively, 
while at 5 AU, k^o = S.ORh, M-.^o = 3.3 M®. 
In a following paper, we will apply these results to 
evaluate whether embryos can attain several earth 
masses while there is adequate residual gas supply 
in the disk for them to acquire their gaseous en- 
velopes and grow into gas giants. 

iv) In the radial velocity surveys, no planet is 
detected in a majority of the target stars. The 
failure for the emergence of any gas giant planets 
does not prevent the embryos to grow after the to- 
tal gas depletion. The eccentricity of the residual 
embryos increases through a post-oligarchic ran- 
dom walk process. As the orbital crossing leads 
to giant impacts, mass growth, and widening sep- 
aration, Tc increases until it is comparable to the 
age of the system. Since Tc is a steeply increas- 
ing function of ko, the separation of embryos is 
unlikely to exceed IORh by much. 

v) However, around stars with known gas giant 
planets, the gas depletion may lead to a sweep- 
ing secular resonance which has the potential to 
shake up the kinematic structure of the "isolated 
embryos" . In Fig. 3b we show that for EMS sys- 
tems which ended up with fco > 10—12, Tc exceeds 
the age of the Solar System. Indeed, the actual 
value of fco is in this range, which accounts for the 
dynamical stability of the Solar System. 

vi) A significant fraction of stars with known 
planets show signs of additional planets. Such 
systems generally have eccentricities much larger 
than those of most planets in the Solar System. 
The emergence of the first-born gas giants in- 
duces the gap formation in their nascent disks and 
the accumulation of planetesimals exterior to the 
outer edge of the gap (Bryden et al. 1999). This 
process promotes the formation of multiple-planet 
systems. In contrast to the embryos, the spac- 
ing between the gas giants may be regulated by 
various migration processes and their masses are 
determining by the disks' thickness-to-radius ra- 
tio. 

Modest ranges of fco and /i values are antic- 
ipated when a system with giant planets forms. 
Gas giants emerging too closely (fco < 5) will un- 
dergo orbital crossing (Fig. 3b), close encounters, 
and cohesive collisions. Gas giants formed with 
fi ~ 10""^ and fco ~ 5.5 have Tc ^ Tdcp whereas 



those with fco ~ 6 have Tc ^ 1 Gyr. The discussion 
under item iii) suggests that close encounters and 
mergers may occur among these gas giant plan- 
ets, which may provide a mechanism for generat- 
ing the large observed eccentricities. We expect 
a considerable dispersion in diffusion rate and the 
asymptotic eccentricities of these systems, because 
gap formation may reduce the efficiency of eccen- 
tricity damping by the planet-disk tidal interac- 
tion. Close encounters between planets with rela- 
tive large masses fi ~ 10"'^ can also lead to non- 
linear effects such as changes of semi-major axis. 
For gas giants formed with fco > 6, neither tidal 
damping nor mutual perturbations of planets are 
effective and they are likely to retain their original 
low-eccentricity orbits. 

vii) We speculate that the large observed ec- 
centricities among the extra solar planets may be 
due to scattering between multiple planets. In §3, 
we show that the asymptotic eccentricities of the 
planets have a Rayleigh distribution, similar to the 
case of planetesimal growth (Ida & Makino 1992, 
Palmer et al. 1993, Lissauer & Stewart 1993). In 
Fig. 13, the eccentricity distribution of the ob- 
served extra solar planets is fit by a Rayleigh dis- 
tribution. The close agreement provides evidence 
that the eccentricity of extra solar planets may be 
excited by the inter-planetary scattering. 

We thank the anonymous referee for valu- 
able suggestions, and Dr. S. Aarseth for im- 
proving the manuscript. This work is supported 
by NSFC(10233020,10778603), NCET (04-0468), 
NASA (NAGS5-11779, NNG04G-191G, NNG06- 
GH45G), JPL (1270927), NSF(AST-0507424, 
PHY99-0794). 



•^We notice after we finished the manuscript that, a similar 
conclusion is also obtained in a recent work by Mario & 
Scott (2007). 
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